1 


Proceedings  -  23rd  Annual  Conference  -  IEEE/EMBS  Oct.25-28,  2001,  Istanbul,  TURKEY 

Tracking  of  nonstationary  EEG  with  Kalman 

smoother  approach 

M.P.  Tarvainen1,2,  P.O.  Ranta-aho1,2  and  P.A.  Karjalainen1 
■'■Department  of  Applied  Physics,  University  of  Kuopio,  Kuopio,  Finland 
2Department  of  Clinical  Neurophysiology,  Kuopio  University  Hospital,  Kuopio,  Finland 


Abstract —  An  adaptive  autoregressive  moving  average 
(ARMA)  modelling  of  nonstationary  EEG  by  means  of 
Kalman  smoother  is  presented.  The  main  advantage  of  the 
Kalman  smoother  approach  compared  to  other  adaptive  al¬ 
gorithms  such  as  LMS  or  RLS  is  that  the  tracking  lag  can 
be  avoided.  This  advantage  is  clearly  presented  with  simula¬ 
tions.  Kalman  smoother  is  also  applied  to  tracking  of  alpha 
band  characteristics  of  real  EEG  during  an  eyes  open/closed 
test.  The  observed  tracking  ability  of  Kalman  smoother, 
compared  to  other  methods  considered,  seemed  to  be  bet¬ 
ter. 

Keywords — Nonstationary  EEG,  Kalman  smoother,  adap¬ 
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I.  Introduction 

Electroencephalogram  (EEG)  analysis  is  a  useful  tool  for 
studying  the  functional  states  of  the  brain  and  for  diagnos¬ 
ing  certain  neurophysiological  states  and  disorders.  In  the 
analysis  of  nonstationary  EEG  the  interest  is  often  to  es¬ 
timate  the  time-varying  spectral  properties  of  the  signal. 
A  traditional  approach  to  this  is  the  spectrogram  method, 
which  is  based  on  Fourier  transformation.  Disadvantages 
of  this  method  are  the  implicit  assumption  of  stationarity 
within  each  segment  and  the  rather  poor  time/frequency 
resolution.  A  better  approach  is  to  use  parametric  spectral 
analysis  methods  based  on  e.g.  time-varying  autoregressive 
moving  average  (ARMA)  modelling.  The  time-varying  pa¬ 
rameter  estimation  problem  can  be  solved  with  adaptive 
algorithms  such  as  least  mean  square  (LMS)  or  recursive 
least  squares  (RLS).  These  algorithms  can  be  derived  from 
the  Kalman  filter  equations  [1],  [2]. 

In  this  paper  we  use  the  Kalman  smoother  algorithm 
in  tracking  of  nonstationary  properties  of  EEG.  Kalman 
smoother  is  compared  to  LMS  and  RLS  algorithms  in 
tracking  of  alpha  band  characteristics  of  EEG  measured 
during  an  eyes  open/closed  test.  The  Kalman  smoother 
approach  is  also  applied  to  the  detection  of  alpha  waves 
of  EEG.  The  main  advantage  of  the  Kalman  smoother  al¬ 
gorithm  compared  to  other  adaptive  algorithms  is  the  fact 
that  the  tracking  lag  can  be  avoided.  This  is  demonstrated 
with  simulations.  Kalman  filter  has  been  previously  used 
in  EEG  analysis  in  e.g.  [3],  [4],  [5]. 

II.  Methods 

If  the  signal  to  be  modelled  is  nonstationary  it  cannot 
be  modelled  as  an  output  of  a  time-invariant  system.  It  is 
natural  in  this  case  to  assume  that  the  system  has  time- 
varying  parameters. 


A.  Time-varying  linear  regression 

Here  we  use  a  time-varying  autoregressive  moving  aver¬ 
age  model  for  the  signal.  The  time-varying  ARMA(p,g) 
model  can  be  written  in  the  form 

p  q 

z(t)  =  -  Y  -  j)  +  Y  bk(t)e(t  ~k)  +  e(t)  (1) 

j— 1  k=l 

where  aj  (t)  and  bk  (t)  are  the  time- varying  ARMA  parame¬ 
ters  and  e(t)  is  the  driving  white  noise  process.  By  denoting 

dt  =  (-ai(i),...  ,-ap{t),bi{t),...  ,bq(t))T  (2) 

(pt  =  {z(t-  1),...  ,z(t-p),e(t-  1),...  ,e(t-q))T( 3) 

the  model  can  be  written  in  the  form 

zt  =  (ft  dt  +  et  (4) 

where  zt  =  z(t)  and  et  =  e(f).  This  is  clearly  a  linear 
observation  model,  with  ipf  being  the  observation  matrix 
and  et  being  the  observation  error.  A  typical  description 
for  the  parameter  variation  when  no  a  priori  information 
is  available,  is  the  random  walk  model  [6].  Thus  for  the 
parameters  9t  we  write  a  state  equation  of  the  form 

9t+ i  =  0t  +  Wt  (5) 

where  Wt  is  a  noise  process.  Equations  (4)  and  (5)  form 
a  specific  form  of  the  general  state  space  equations,  with 
the  input  process  wt-  Now  the  problem  is  to  estimate  the 
time- varying  parameters  9tl  according  to  the  state  space 
model. 

B.  Kalman  filter 

The  Kalman  filtering  problem  is  to  find  the  minimum 
mean  square  estimator  9t.  for  state  9t  given  the  observa¬ 
tions  z\, . . .  ,zt-  This  has  been  shown  to  be  equal  to  the 
conditional  expectation  value 

9t  =  E{9t\zi,...  ,zt}  (6) 

We  assume  here  the  state  and  measurement  noises  Wt  and 
et  to  be  uncorrelated,  zero  mean,  random  processes  with 
covariance  matrices  CWt  =  cr2/  and  Cet  =  cr2/,  so  that 
the  individual  parameter  evolutions  are  assumed  to  be  in¬ 
dependent.  The  initial  state  9q  is  assumed  to  be  uncor¬ 
related  with  et  and  Wt  with  finite  variance.  The  Kalman 
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where  9t\t-i  is  the  mean  square  estimator  for  state  Ot  given 
the  observations  z±, . . .  ,  Zt~  1,  0t  is  the  state  estimation  er¬ 
ror  6t  =  9t  —  9t  and  Kt  is  the  Kalman  gain  matrix.  The 
adaptation  of  the  filter  is  primarily  affected  by  CWt.  Note 
that  the  unknown  process  et  is  here  estimated  with  the 
prediction  error  process  e*  in  every  step  of  the  iteration. 

C.  Fixed-interval  smoother 

The  fixed-interval  smoothing  problem  is  to  determine  es¬ 
timates 


(b) 


9t\T  —  E  {9t\zi, . . 

. .  , zT} 

(13) 

for  fixed  T  and 
solution  for  this 

for  all  t  in  the  interval  1  <  t  <  T. 
can  be  written  in  the  form  [7] 

The 

@t—l\T 

=  9t-i+At-i 

1 

T 

(14) 

=  cs  C-1 

“t-i  0t  |t_1 

(15) 

where  At_i  includes  the  error  covariances  stored  in  the  for¬ 
ward  run  of  Kalman  filter.  Also  the  state  estimates  9t  and 
9t u  i  need  to  be  stored.  The  smoothed  estimates  9t_ are 
then  obtained  by  running  the  stored  estimates  backwards 
in  time  by  taking  t  =  T,  T  —  1, . . .  ,2.  The  initialization  is 
evidently  with  the  filtered  estimate. 


Fig.  1.  AR(2)  process  estimation  with  RLS  and  Kalman  smoother 

algorithms.  The  root  evolution  and  the  realization  are  presented 
in  block  (a).  Both  algorithms  were  optimized  (b).  Optimal  value 
for  the  state  noise  covariance  coefficient  of  the  Kalman  smoother 
was  (J^,  =  0.001  and  the  forgetting  factor  of  RLS  was  A  =  0.935. 
The  estimates  of  the  modulus  and  phase  angle  of  the  root  are 
shown  in  block  (c).  The  true  values  (black),  Kalman  smoother 
estimates  (red)  and  optimal  RLS  estimates  (blue).  The  smoother 
RLS  estimates  (green)  were  calculated  by  using  A  =  0.98. 

III.  Results 

In  order  to  evaluate  the  performance  of  the  Kalman 
smoother  algorithm  we  conduct  two  simulations,  where 
Kalman  smoother  is  compared  to  the  popular  forgetting 
factor  RLS  algorithm.  Finally  the  Kalman  smoother  is  ap¬ 
plied  to  time- varying  spectrum  estimation  of  real  EEG  and 
for  alpha  wave  detection. 


D.  Spectral  estimation 


A.  Simulations 


Once  the  time-varying  coefficients  of  the  ARMA(p,g) 
model  (1)  are  solved  the  time- varying  power  spectral  den¬ 
sity  (PSD)  estimation  can  be  obtained  in  the  terms  of  the 
estimated  coefficients 


PtW 


Kt) 


|i  +  EtAWe—T 


(16) 


where  cr2(t)  is  the  prediction  error  variance.  After  the 
adaptive  algorithm,  used  to  estimate  the  time- varying 
ARMA  parameters,  converges  power  spectrum  can  be  cal¬ 
culated  for  each  time  instant.  However  in  real  applications 
it  is  not  usually  necessary  to  calculate  PSD  estimate  for 
each  time  instant,  but  rather  only  once  for  a  certain  time 
interval.  The  average  spectrum  within  certain  interval  can 
be  calculated  from  averaged  ARMA  parameters  of  that  in¬ 
terval. 


In  the  first  simulation  a  time-varying  signal  was  gener¬ 
ated  as  an  AR(2)  process.  The  root  evolution  and  a  typ¬ 
ical  realization  are  presented  in  Fig.  1  (a).  The  modulus 
and  phase  angle  of  the  root  were  estimated  with  Kalman 
smoother  and  RLS  algorithms.  Parameters  controlling  the 
adaptation  were  optimized  in  both  algorithms  to  obtain  the 
minimum  error  in  AR  coefficient  estimation.  The  estima¬ 
tion  errors  as  a  function  of  adaptation  parameters  for  both 
algorithms  are  presented  in  Fig.  1  (b).  The  estimates  are 
shown  in  Fig.  1  (c). 

Two  RLS  estimates  were  calculated  to  demonstrate  the 
effect  of  the  forgetting  factor  on  the  estimates.  RLS  esti¬ 
mates  with  the  optimal  value  for  the  forgetting  factor  have 
only  a  small  tracking  lag  but  the  estimates  are  far  more 
unstable  compared  to  the  Kalman  smoother  estimates.  By 
increasing  A  RLS  estimates  become  more  stable  but  the 
tracking  lag  increases  at  the  same  time.  This  simulation 
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Fig.  3.  Time-varying  spectral  analysis  of  ERD/ERS  test  of  alpha 
waves  of  EEG.  The  measured  EEG  from  channel  02  is  shown 
on  the  topmost  axis.  The  time  window  used  in  the  spectrogram 
was  2  seconds.  The  step  size  of  LMS  algorithm  was  //  =  0.0002 
while  the  forgetting  factor  of  RLS  was  A  =  0.95.  The  state  noise 
covariance  coefficient  of  the  Kalman  filter  was  <r?.  =  0.0003. 


Fig.  2.  A  realistic  simulation  of  EEG  transition  as  an  AR(5)  pro¬ 
cess.  (a)  The  roots  and  the  corresponding  spectra  before  (blue) 
and  after  (red)  the  transition,  (b)  A  typical  realization  of  the 
process.  Averaged  estimates  over  100  realizations  of  the  modu¬ 
lus  and  phase  angle  of  the  root  corresponding  to  alpha  activity 
are  presented  in  (c) ,  where  true  values  (black) ,  Kalman  smoother 
estimates  (red)  and  RLS  estimates  (blue/green)  are  shown.  The 
state  noise  covariance  coefficient  of  the  Kalman  smoother  was 
crj,  =  8  ■  10— 5  and  the  forgetting  factors  of  RLS  were  Ai  =  0.98 
(blue)  and  A2  =0.9  (green). 

shows  clearly  the  advantages  of  the  Kalman  smoother  com¬ 
pared  to  the  RLS  algorithm.  However  not  much  can  be  said 
about  the  performance  of  the  Kalman  smoother  in  track¬ 
ing  of  nonstationary  EEG  based  on  this  simple  simulation. 
Hence  we  aim  to  a  more  realistic  simulation  of  EEG. 

In  many  cases  we  are  interested  in  tracking  of  nar¬ 
row  band  characteristics  of  the  EEG  signal.  One  such 
case  is  the  event  related  desynchronization/synchronization 
(ERD/ERS)  of  alpha  waves.  The  occipital  EEG  recorded 
while  patient  having  eyes  closed  shows  high  intensity  in 
the  alpha  band  (7-13  Hz).  With  the  opening  of  the  eyes 
this  intensity  decreases  or  even  vanishes.  It  can  be  as¬ 
sumed  that  EEG  exhibits  a  transition  from  a  stationary 
state  to  another.  Such  a  transition  was  here  simulated 
as  an  AR(5)  process.  The  roots  of  the  system  for  both 
stationary  states  (obtained  from  real  EEG  measurements) 
and  the  corresponding  power  spectrums  are  presented  in 
Fig.  2  (a),  where  the  strong  peak  around  0.387T  is  due  to 
the  synchronization  of  the  alpha  waves  of  EEG. 

In  order  to  make  the  simulation  more  realistic  abrupt 
transitions  of  AR  coefficients  were  smoothed  as  described 


in  [8].  A  typical  realization  of  the  simulated  AR(5)  process 
is  presented  in  Fig.  2  (b).  Results  of  tracking  the  alpha 
band  characteristics  are  presented  in  Fig.  2  (c),  where  av¬ 
eraged  estimates  over  100  realizations  of  the  phase  angle 
and  the  magnitude  of  the  root  corresponding  to  alpha  ac¬ 
tivity  are  presented.  In  order  to  obtain  as  smooth  estimates 
with  RLS  as  is  obtained  with  Kalman  smoother  the  forget¬ 
ting  factor  A  must  be  quite  small.  However  this  leads  to 
substantial  tracking  lag.  With  larger  values  of  A  the  track¬ 
ing  lag  can  be  attenuated,  but  estimates  become  now  more 
unstable. 

B.  ERD/ERS  of  alpha  waves  of  EEG 

The  eyes  open/closed  test  is  a  typical  application  of  test¬ 
ing  the  desynchronization/synchronization  of  alpha  waves 
of  EEG.  One  such  transition  from  desynchronized  state  to 
synchronized  state  is  presented  in  Fig.  3.  The  performance 
of  the  Kalman  smoother  in  tracking  of  alpha  band  charac¬ 
teristics  is  compared  to  other  commonly  used  adaptive  al¬ 
gorithms  such  as  RLS  and  LMS  and  also  to  the  traditional 
spectrogram  method.  An  ARMA(6,2)  model  was  used  in 
all  adaptive  algorithms.  The  length  of  the  time-window 
used  in  spectrogram  was  2  seconds,  which  is  long  enough 
when  considering  the  frequencies  of  the  alpha  band  (7-13 
Hz).  The  step  size  of  the  LMS  algorithm  was  /r  =  0.0002 
and  the  forgetting  factor  of  RLS  was  chosen  to  be  A  =  0.95 
resulting  in  quite  stable  estimates  and  still  rather  fast  adap¬ 
tivity.  The  state  noise  covariance  coefficient  of  the  Kalman 
smoother  was  crj  =  0.0003.  The  tracking  speed  of  the 
Kalman  smoother  seems  to  be  fastest  and  an  interesting 
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The  Kalman  smoother  algorithm  was  applied  to  track¬ 
ing  of  nonstationary  EEG.  The  performance  of  Kalman 
smoother  in  tracking  of  alpha  band  characteristics  seemed 
to  be  most  reliable  compared  to  LMS  and  RLS  algorithms. 
Kalman  smoother  was  also  applied  to  the  detection  of  alpha 
waves  of  EEG  with  success.  Also  two  simulations  were  con¬ 
ducted  showing  clearly  the  main  advantages  (smooth  esti¬ 
mates  without  tracking  lag)  of  Kalman  smoother  compared 
to  other  adaptive  algorithms.  The  implementation  and  us¬ 
ability  of  the  Kalman  smoother  approach  are  straightfor¬ 
ward.  The  adaptation  rate  is  adjusted  simply  by  setting 
the  state  covariance  coefficient 

One  problem  in  modelling  the  data  with  adaptive  al¬ 
gorithms  is  the  selection  of  the  model  order.  For  time- 
invariant  systems  there  exist  various  criteria  for  selecting 
the  model  order  [9] .  All  these  criteria  are  based  on  the  com¬ 
promise  between  model  fit  and  model  complexity.  Also  in 
the  time-varying  case  there  exist  some  criteria  for  selecting 
the  model  order.  For  example  in  [10]  the  use  of  Akaike’s  in¬ 
formation  criterion  (AIC)  was  justified  in  the  time- varying 
case  under  certain  conditions.  However  in  the  case  of  track¬ 
ing  alpha  rhythm  of  EEG  the  ARMA  model  of  order  p  =  6 
and  q  =  2  seems  to  be  suitable.  The  same  model  order  was 
also  used  in  [11],  [12]. 


Frequency  (Hz) 

(c) 

Fig.  4.  Kalman  smoother  applied  to  alpha  rhythm  detection,  (a) 
An  EEG  sample  of  15  seconds  from  channel  02  measured  from 
subject  having  eyes  closed  and  the  corresponding  time-varying 
PSD.  (b)  Detection  is  based  on  thresholding  the  power  integral 
over  the  alpha  band  (7-13  Hz)  with  a  threshold  of  10  /xV2/ Hz. 
Block  (c)  presents  PSD  estimates  (calculated  with  traditional 
FFT  based  method)  for  the  signals  obtained  by  concatenating 
the  EEG  epochs  where  alpha  activity  was  detected  (red)  or  not 
detected  (blue). 

gap  in  alpha  rhythm  is  observed  after  9  seconds.  The  con¬ 
tents  of  this  kind  of  gaps  is  considered  in  Fig.  4. 

C.  Detection  of  alpha  rhythm  of  EEG 

The  aim  of  automatic  EEG  analysis  is  often  the  detec¬ 
tion  of  certain  waveforms.  Hence  the  performance  of  the 
Kalman  smoother  on  detection  of  alpha  waves  of  EEG  is 
considered  here.  Fig.  4  (a)  presents  a  time- varying  spec¬ 
trum  for  an  EEG  sample  of  15  seconds  measured  from 
healthy  subject  having  eyes  closed.  Alpha  wave  detection 
was  obtained  by  thresholding  the  power  integral  over  the  al¬ 
pha  band  (7-13  Hz).  The  threshold  was  set  to  10  /iV2/Hz. 
The  power  integral  curve  was  filtered  with  a  median  filter 
of  third  order  to  obtain  smoother  detection.  The  perfor¬ 
mance  of  the  alpha  detector  was  verified  by  concatenating 
the  EEG  epochs  where  alpha  waves  were  detected  and  those 
were  no  detection  was  made.  The  PSD  estimates,  calcu¬ 
lated  with  a  traditional  FFT  based  periodogram  method, 
for  these  concatenated  signals  are  presented  in  Fig.  4  (c) 
verifying  the  absence  of  alpha  rhythm  in  the  lower  concate¬ 
nated  signal. 
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